Glassiness, Rigidity and Jamming of Prictionless Soft Core Disks 
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The jamming of bi-disperse soft core disks is considered, using a variety of different protocols to 
produce the jammed state. In agreement with other works, we find that cooling and compression 
can lead to a broad range of jamming packing fractions <f>j, depending on cooling rate and initial 
configuration; the larger the degree of big particle clustering in the initial configuration, the larger 
will be the value of cj>j. In contrast, we find that shearing disrupts particle clustering, leading to a 
much narrower range of <j)j as the shear strain rate varies. In the limit of vanishingly small shear 
strain rate, we find a unique non-trivial value for the jamming density that is independent of the 
initial system configuration. We conclude that shear driven jamming is a unique and well defined 
critical point in the space of shear driven steady states. We clarify the relation between glassy 
behavior, rigidity and jamming in such systems and relate our results to recent experiments. 
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I. INTRODUCTION 

An athermal system of hard or soft core interacting 
particles, for which thermal fluctuations are negligible 
(i.e. T = 0), is found to undergo a jamming transition 
from a liquid-like state to a rigid but disordered solid 
as the packing fraction <f> increases above a critical value 
4>j pp. Other physical systems similarly seem to show a 
transition from a liquid to a rigid but disordered state, 
as a function of some physical control parameter: foams 
change from flowing liquid to elastic solid once the ap- 
plied shear stress a falls below a critical shear yield stress 
cry; supercooled liquids appear to freeze into a frozen 
glass as the temperature T is lowered. 

These observations led Liu, Nagel and co-workers [2J [3J 
to propose that these phenomena might be unified in 
terms of a three dimensional jamming phase diagram 
with the axes of packing fraction <f>, applied shear stress a, 
and temperature T . A surface in this three dimensional 
phase space separates jammed (i.e. rigid but disordered) 
from unjammed (i.e. liquid-like) states. We sketch this 
jamming phase diagram in Fig. [Tl following Ref. [2j. As 
originally proposed [2J [3J, this jamming surface repre- 
sented points in phase space where the relaxation time r 
of the system reached some experimentally defined large 
time scale. The intersection of this surface with the equi- 
librium ((f), T) plane at a = then describes the exper- 
imentally observed glass transition T° xp (0), where the 
viscosity of the liquid becomes immeasurably large upon 
cooling. 

For discussing possible critical behavior, it is of inter- 
est theoretically to consider the critical jamming surface 
that would correspond to a diverging time scale r — > oo. 
In this case the line T g ((f) in Fig. HI would be a true 
equilibrium glass transition. By saying equilibrium glass 
transition, we mean that T g (<f>) would locate a true sin- 
gularity in the free energy, independent of the particular 
dynamics of the system. Subsequently, when we refer to 



equilibrium glass transition, this will be what we mean. 
Such a critical jamming surface would intersect the <f> 
axis at T = a = at a well defined point </>o, such that 
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O'Hern et al. [2J conjectured 



this point to be identical to the T — jamming tran- 
sition of athermal particles, i.e. </>o = <pj, and denoted 
it "point J". Moreover, they conjectured that point J 
may act like a critical point and "that it may control the 
region around it and thereby govern the nature of the 
entire jamming surface in the phase diagram" (2j. Were 
this conjecture correct, then properties of the equilibrium 
glass transition would be intimately related to properties 
of athermal jamming. In the following, we will denote 
"point J" as the athermal jamming point, within some 
well defined physical protocol, where the packing fraction 
equals </>,/ and T = a = 0. 




point J 



1/0 surface below which 
states are jammed 

FIG. 1: Proposed jamming phase diagram, adapted from 
Ref. [2]. Note that the axis along which packing fraction 
varies is l/<f> rather than <j>. 



It must be noted that in models of simple liquids, 
such as those considered in this work, the existence of 
an equilibrium glass transition at a finite temperature, 
Tg(tfi) > 0, remains a much debated question. Works 



have suggested that the critical jamming surface, where 
t — > oo, may collapse entirely into the T — > plane. See 
Ref. 4 for a recent review. However, even in this case, 
one might expect at a = a line of zero temperature glass 
transitions, that would terminate at some specific lowest 
packing fraction cf>o, with viscosity diverging as T — > for 
all <j) > <Pq. If such a 0o was identical to the <fij of point J, 
one would again have a connection between equilibrium 
glassy behavior (albeit a T g — > glass transition) and 
athermal jamming. 

Although the above conjecture is appealing, several re- 
cent works have suggested that the actual situation may 
be more complicated. Some simulations have suggested 
that an equilibrium glass transition may be controlled by 
a different critical point, sometimes referred to as "point 
G", that is distinct from the athermal jamming point J. 
Equilibrium simulations of hard spheres (where tempera- 
ture T plays no role and density (f> is the only parameter) 
in three dimensions (3D) by Brambilla et al. 5 claim a 
glass transition at a 4>q that is distinctly lower than the 
typical values of 4>j obtained from compression; the re- 
duced pressure p/T remains finite at <fio, in contrast to the 
athermal jamming of hard spheres where p/T is expected 
to diverge. Equilibrium simulations of soft spheres in 3D 
by Berthier and Witten [6l[7] show a scaling in the (</>, T) 
plane that similarly suggests aT->0 glass transition at 
a 0o lower than the athermal jamming <pj. Similar re- 
sults have been suggested by Xu et al. (SJ. Other works 
by Donev et al. [9] argue that in a properly equilibrated 
hard sphere system there is no glass transition, unless 
some constraint is imposed to prevent crystallization. 

At the same time, recent works have illustrated that 
the precise value <f>j of the athermal point J is not unique, 
even within a given specified model, but can depend upon 
the particular protocol used to prepare the jammed state. 
Donev et al. [9] show, for both frictionless mono-disperse 
spheres in 3D and bi-disperse disks in two dimensions 
(2D), that for compression driven jamming, (j)j depends 
on the rate of compression. Chaudhuri et al. [TU] show 
that when compressing configurations equilibrated at an 
initial </>i n ;t, bi-disperse frictionless spheres in 3D jam at 
a (jjj that increases with </>i n it; hence the value of 4>j can 
depend not only on the compression rate, but on the 
particular initial configuration from which one starts the 
compression. Recent theoretical works [TTJ Q~2] on mean- 
field hard sphere models have found similar results: a 
continuous range of athermal jamming densities at infi- 
nite p/T, with a cj)j that varies according to compression 
or cooling rate, as well as a distinct equilibrium glass 
transition at a finite p/T. 

In this work we consider a variety of jamming proto- 
cols for a 2D system of frictionless bi-disperse soft-core 
disks, focusing on protocols which do not involve com- 
pression. We also carry out equilibrium Monte Carlo 
(MC) simulations of hard disks to look for the onset of 
glassy behavior prior to jamming, i.e. at a (f>o < 4>j, as 
was previously observed in 3D. Our main conclusion is 
that, unlike jamming driven by compression or cooling, 



athermal shear driven jamming results, in the limit of a 
vanishingly small shear strain rate 7 — > 0, in a unique, 
well-defined, non-trivial, value of <f>j that is independent 
of the system's starting configuration. This result follows 
from our observation that shearing breaks up the cluster- 
ing of big particles that can lead to phase separation and 
crystallization under slow compression or cooling. The 
remainder of this paper is organized as follows. In Sec- 
tion II we describe our model of bi-disperse frictionless 
soft core disks, give our precise procedure for determining 
the jamming packing fraction </>j, and describe the dif- 
ferent jamming protocols we will consider. In Section III 
we present our numerical results. In Section IV we relate 
our results to some recent experiments and summarize 
our conclusions. 



II. MODEL 

Our system is a bi-disperse mixture of frictionless disks 
with diameter ratio g?b/<^s = 1-4, as has been used in 
earlier works O [13] . The fraction of bigger particles is 
ib = 1/2. The disks interact via a pairwise soft-core re- 
pulsive contact interaction, that is harmonic in the par- 
ticle overlap, 
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where Vij is the distance between the centers of two par- 
ticles i and j, and dij is the sum of their radii. Length 
is measured in units such that the smaller diameter is 
unity, and energy is measured in units such that e = 1. 
A system of N disks in an area A thus has a packing 
fraction (density), 



<j) = Nir{0.5 2 +0.7 2 )/(2A) 



(2) 



We use a simulation cell of area A = L x L y , with equal 
length and width, L x = L y . 

Starting from a set of physically motivated initial 
states at fixed <j), we use a non-linear conjugate gradient 
method to quench each state to its local energy mini- 
mum (the inherent structures). For soft core particles, 
mechanically stable states exist at values of 4> above the 
jamming cf>j. Such states are characterized by a finite 
interaction energy, pressure, and shear yield stress; all 
these vanish as <f> — > <f>j from above [3j- Energy mini- 
mized states with an energy per particle below a certain 
fixed very small threshold value, E/N < e cut , are there- 
fore regarded as unjammed; otherwise the state is con- 
sidered to be jammed. This criterion for jamming [2] has 
been shown [14] to be essentially the same as "strictly 
jammed" in the classification scheme of Torquato and 
Stillinger [TS]. In this manner we count the fraction /(</>) 
of these energy minimized states which are jammed. 

As <j> increases, /(</)) varies rapidly from zero to unity, 
with /(</>) approaching a sharp step function as the num- 
ber of particles N — > 00. The location of the step then 



determines the jamming packing fraction cf>j for that ini- 
tial set of states. We consider two classes of initial states: 
(i) Equilibration at a finite T, which may be thought of 
as the equilibrium temperature of a glassy system, or as 
an effective temperature of kinetic motion in a granular 
system with uniform mechanical agitation. Quenching 
corresponds to suddenly turning off the agitation and al- 
lowing the system to relax. In the limit T — oo one 
chooses random initial positions. This is the ensemble 
studied by O'Hern et al. [2] and we will denote it as 
"RAND." (ii) Shearing at a constant uniform shear strain 
rate 7. Quenching corresponds to suddenly turning off 
the shear and allowing the system to relax to a mechan- 
ically stable or to an unjammed state. The limit 7^0 
gives quasistatic shearing ("QS"), as studied previously 
by Heussinger and co-workers |16l I17j . 

The specific conjugate gradient algorithm we use to 
energy minimize is the Polak-Ribiere method |18j . We 
stop the energy minimization when one of the following 
conditions is met: (i) the relative decrease in the energy 
AE/E after 50 iterations is smaller than 10~ 10 , or (ii) the 
energy per particle falls below a certain small threshold 
value, E/N < e cut . In the second case, we consider the 
state to be unjammed. We find that the threshold value 
e c .ut = 10 -16 gives a clear separation between the jammed 
and unjammed states up to the largest system size we 
have studied (see appendix A for further details). 



III. RESULTS 

A. Random vs Quasistatic Shearing Ensembles 

In Fig. [2] we plot our results for the RAND and QS en- 
sembles, showing how the jammed fraction /(</>) sharpens 
to a step function as the number of particles N increases. 
For RAND we average over at least 10000 initial configu- 
rations for each value of </>. For QS we average over 10—20 
independent runs, each sheared a total strain 7^4 — 8 
for our biggest size, but 7 = 200 for our smallest size. 
We use Lees-Edwards boundary conditions [19] to model 
uniform shear strain, energy minimizing after each small 
strain step A7. For N > 4096 we use A7 = 10~ 5 , while 
for N < 4096 we use A7 = 10~ 4 . We have explicitly ver- 
ified that for all, but perhaps the biggest size N — 8192, 
our value of A7 is small enough not to influence our re- 
sults (see appendix A for details). We clearly see that 
the two ensembles approach different jamming densities, 
0.842 while <t>J S ~ 0.843 [2D]. 



^RAND 



B. Equilibration at Finite Temperature 

Next we consider initial states equilibrated at a fixed 
temperature T. In Fig. [3k we plot the jammed fraction 
f(<p) resulting from the energy minimized states arising 
from these thermally equilibrated initial states, compar- 
ing RAND with several finite values of T, for N = 256 
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FIG. 2: (color online) Jammed fraction / vs packing fraction 
4> for systems of different numbers of particles TV. (a) and 
(b) are for the RAND and QS ensembles, respectively. Ver- 
tical dashed lines indicate the limiting TV — > 00 value of the 
jamming density cf>j in each case. 



particles. For the three lowest T we also show results 
for N = 512 to illustrate that increasing N continues to 
lead to a sharpening of the transition as seen in Fig. [2] 
To equilibrate at T we do ordinary MC simulations, at 
each step displacing a randomly chosen particle by a ran- 
dom amount and accepting or rejecting the move accord- 
ing to the Metropolis algorithm. N attempted particle 
moves is 1 MC pass. At our lowest T, we use 10 in- 
dependent runs, each with roughly 10 8 MC passes. We 
judge that we have equilibrated when particles have, on 
average, diffused a distance equal to a few particle diam- 
eters. We see that T = 5 x 10~ 3 is essentially equal to the 
T = 00 ensemble RAND, but that as T decreases, </>j(T) 
increases. Our lowest T* = 5 x 10~ 4 gives our largest 
4>j(T*) ~ 0.850 [2TJ . For such high densities, our runs at 
T* are just at the border of equilibrating; we would need 
much longer runs to try and equilibrate at even lower T. 
Similar results have recently been found for continuous 
cooling with different fixed rates [22 . These results are 
in good agreement with recent predictions from a mean- 
field-like hard-sphere model ill . 

Roughly the same range of <pj was found by Donev et 
al. [H] from slow compression of bi-disperse hard disks 
(they use xb — 1/3). Donev et al. argue that an in- 
creased <pj results from an increased order due to the 
clustering of big particles. To check for such clustering 
we have computed the fraction, ng, of big particles which 
have 6 nearest neighbors (as determined by Delaunay tri- 
angularization) that are also big particles. In Fig. [3b we 
plot n 6 vs (j) for the cases of Fig. [3b, as well as QS. We 
see little difference in uq comparing RAND, QS, and the 
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FIG. 3: (color online) (a) Jammed fraction / vs packing frac- 
tion <f> comparing RAND with ensembles quenched after ther- 
mal equilibration at a fixed temperature T. Open symbols 
(solid lines) are for TV = 256 particles, while closed sym- 
bols (dashed lines) are for N = 512 particles, (b) Clustering 
strength parameter n& vs cj> for the cases shown in (a) as well 
as for QS. Representative error bars are shown at select data 
points at the lowest T's. 



highest T, however n§ systematically increases as T de- 
creases. The increasing fluctuations in n 6 as <fi varies 
at low T reflect the increasing difficulty to equilibrate. 
Donev et al. have argued that, given sufficiently long 
equilibration, even higher values of 4>j might be achieved, 
up to the maximum max ~ 0.91 of fully phase separated 
lattices of big and small particles. We expect a similar 
situation in our system, if we could equilibrate at even 
lower T. 



C. Shear Driven Steady States 

Next we consider shearing the system. For our simula- 
tions at a finite shear strain rate 7 we use Durian's [23] 
foam dynamics: overdamped motion with viscous damp- 
ing to the local average shear flow velocity. For shear 
flow in the x direction we have, 



them and count the resulting fraction that are jammed. 
In Fig. Rlwe show results for the jammed fraction f((f>) 
for TV = 256 particles. 

We stress that the energy minimization step, repre- 
senting a sudden switch off of the applied strain rate 7, is 
crucial to this calculation. Were we to sample the steady 
state distribution at fixed 7 directly, all states would have 
a finite energy yet all states are flowing; our criterion 
for computing the jammed fraction is only applicable to 
mechanically stable states at rest. Thus the 4>j(j) de- 
termined by the present procedure should not be viewed 
as representing a jamming transition for driven steady 
states at finite 7 (there is none, since all such states by 
definition are flowing); rather it represents the jamming 
point resulting from a particular dynamic protocol for 
creating statically jammed states, i.e. relaxation to rest 
from initial states driven at a finite shear rate 7. 

Our fastest shear rate 7 = 1CP 3 gives results equal 
to the random initial positions of RAND. From this we 
can infer that, in a rapidly sheared system, the soft core 
interactions are playing little role in ordering the parti- 
cles, and the system is thus passing through effectively 
random configurations. Our slowest shear rate 7 = 10~ 8 
is clearly converging to the QS limit. Thus quasistatic 
shearing, in which the system is always instantaneously 
relaxing into its nearest local energy minimum as it is 
slowly sheared, is the appropriate 7 — > limit of the 
overdamped dynamics of Eq. (pi). From our results in 

Figs. [2p and 12] we thus conclude that there is indeed a 

— • — • • OS 

well defined jamming density (pj in the 7 — > 0, N — > 00, 

limit. 
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FIG. 4: (color online) Jammed fraction / vs packing fraction 
(f> for a system of TV = 256 particles, comparing RAND and 
QS ensembles with ensembles quenched from a fixed finite 
shear strain rate 7. 
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D. Hard Disk Equilibrium Monte Carlo 



where the last term yci x is the average shear flow veloc- 
tiy. Lees-Edwards boundary conditions are used, and we 
choose units of time such that C = 1. We run the sim- 
ulations up to a certain total strain 7 = jt (7 = 33 for 
our smallest 7), and then sampling configurations uni- 
formly from the resulting shear flow, we energy minimize 



Next we consider the dynamic behavior at low densities 
4> < </>j AND . Following Brambilla et al. jS], we simulate 
the diffusion of N = 1024 hard core disks using local MC 
moves in which a randomly selected particle is displaced 
a random amount within a box of length O.lds about 
its center; a move is accepted only if the non-overlap 



hard disk constraint is obeyed. N such attempted moves 
corresponds to one MC pass, which we equate to one 
unit of time. Measuring the self-part of the intermediate 
scattering function [5], 



F s ( g ,^i£>^«- 



r<(0)] 



(4) 



with q = (2w/1.2ds)x, we define the relaxation time r 
by F s (q, r) = 1/e. In Fig. [5h we show results for r vs 0. 
At low (f> < 0.76 equilibration is relatively straight for- 
ward. At larger <fr, we use the following procedure to try 
and stay on the metastable glassy part of the equation 
of state: starting from the ending configuration of the 
previous value of 4>, we compress the system an amount 
Acf> — 0.005, and then simulate for a time of roughly 
100r before increasing cf> again. At our lowest 4> this cor- 
responds to 3 x 10 5 MC passes; for our highest <f> this is 
5 x 10 8 MC passes. We leave aside the question whether 
t(4>) is truly diverging at an ideal glass transition (f>o, 
as suggested by Berthier and Witten [7] , or whether 
the growth in r is a kinetic effect of falling out of equi- 
librium, as argued by Donev et al. [5]. Here we just 
note that r clearly grows many orders of magnitude by 
the time one reaches O ~ 0.80 < <^ AND ~ 0.842, thus 
leading to glassy behavior before the onset of our lowest 
jamming density. In Fig. [5b we show the corresponding 
cluster strength ne, which we see increases rapidly with 
increasing <\>. 

In their work, Chaudhuri et al. [10] observed that when 
they equilibrated their system first as hard spheres at 
some initial </>i n itj and then slowly compressed them to 
reach jamming, the 4>j that resulted increased with in- 
creasing </>i n it. From Fig.[5]D we see that, for equilibrated 
systems, the clustering strength n$ increases rapidly with 
increasing cf>. The results of Chaudhuri et al. are thus 
consistent [24] with the assertion by Donev et al. [9] 
that initial configurations with greater clustering result 
in jamming at larger </>j, when compressed. 
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FIG. 5: (a) Relaxation time r and (b) clustering strength 
parameter n% vs <fi for N = 1024 diffusing hard core particles. 
Representative error bars are shown. 



The data of Fig. [5] does not represent true equilibrium 
at the largest values of <fi shown. We have found that 
we are able to more properly equilibrate the system if we 



include non-local swaps between big and small particles 
in our MC moves. Such moves are, of course, unphysical 
when modeling a continuous dynamics of the particles, 
but they are perfectly acceptable for sampling true equi- 
librium. With such moves we find evidence for a tran- 
sition near 4> ~ 0.78 to a phase separated coexistence 
between a liquid mixture of big and small particles and a 
solid of big particles, just as was predicted by Donev et 
al. Pj. In such true equilibrium states, n^ becomes even 
larger than found in Fig. [5b. 



E. Compression vs Shear Driven Jamming 

To illustrate our above results on jamming, we next 
consider the following numerical "experiment" . Since the 
largest values of n§ in Fig. [SJa are slightly larger than 
found in Fig. [3j) from cooling soft disks, we expect that 
compression of configurations equilibrated at densities 
4> ~ 0.80 should result in relatively high jamming den- 
sities. We therefore take one configuration at (f> = 0.80, 
sampled from the states that produced the data of Fig. [5] 
we denote this as configuration "A". We take a second 
configuration "B", obtained also at <j> = 0.80, but after 
doing MC with particle swaps so as to achieve a better 
equilibration of the system and a higher degree of parti- 
cle clustering. Configuration A has n$ = 0.037 while B 
has n6 = 0.168. Both have N = 1024 particles. We show 
these initial configurations A and B in Fig. [6J 




(a) configuration A 



0.80 (b) configuration B 



FIG. 6: Initial unjammed configurations A and B at <j> — 0.80. 
Configuration B is seen to have a greater degree of big (black) 
particle clustering than A. N = 1024. 



We then uniformly compress both configurations A and 
B in steps of A<p = 10~ 4 , relaxing the system to its local 
energy minimum after each compression step. In Fig. [7k 
we plot the resulting energy per particle E/N vs 4>- We 
see that A jams at the relatively high value of 0j ~ 
0.8534, while B jams at the even higher </>j ~ 0.8559. 
We then return to these configurations as they were at 
4> = 0.85. Because 0.85 is below the jamming density 
of either system, these are unjammed, stress-free, states. 
We now quasistatically shear these configurations using a 
strain step A7 = 10~ 4 . Our results are shown in Fig.[7p. 



We see that after relatively small strains of 7 = 0.05 for 
A, and 7 = 0.077 for B, both systems jam. This is as 
expected since <fi = 0.85 > <fij . In Fig. [7c we plot the 
instantaneous value of n§ for these two configurations, as 
a function of total shear strain 7 at <fi — 0.85, showing 
results out to a much larger total strain 7 = 35 than in 
Fig. [Tfj. We see that after a certain amount of shearing, 
Uq for both A and B drop down to the values typical of 
the QS ensemble (see Fig. [3]b). We thus conclude that 
shearing breaks up the particle clustering that can lead 
to high jamming densities under slow compression, and 
that this effect is responsible for the well defined jamming 
density cpj in the limit of vanishingly small shear rate, 
7->0. 
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FIG. 7: (color online) (a) Energy E/N vs <fi for the two config- 
urations A and B of Fig. [6] undergoing uniform compression. 

(b) E/N vs strain 7 for the same configurations undergoing 
uniform shear starting from stress-free states at <fi = 0.85. 

(c) Clustering strength n§ vs 7 as the two configurations are 
sheared at d> = 0.85. N = 1024. 



To further illustrate this point, in Fig. [5] we plot the 
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the strain averaged energy per particle (E/N) 7 , and the 
strain averaged jammed fraction (/) 7 vs total shear strain 
7, as we quasistatically shear three different initial states: 
the moderately clustered state A and the highly clustered 
state B, as in Figs. [6] and [7J and a state C that starts with 
particles in random positions. 

In Fig. [8k,b,c, we show results at the density <fi = 
0.85 > 4>j (the same density as in Fig. r7p). In 
Figs. [8H,e,f we show results at the lower density <fi — 

— OS 

0.84 < cpj .In both cases we see that as 7 increases, the 
strain averaged quantities for the different initial con- 
figurations approach a common steady value. For the 
moderately clustered initial state A, and the random ini- 
tial state C, this happens after a relatively short strain; 
for the highly clustered initial state B, it takes consider- 
ably longer to lose memory of the initial state. Note, for 
<fi = 0.84, / sw 0.05 and so the system is rarely jammed; as 
a finite energy comes only from jammed configurations, 
the rarity of jammed configurations at the low <fi = 0.84 
is the reason for the larger fluctuations observed in the 
curves for (E/N)^ shown in Fig. [8k 

These observations illustrate two important points: (i) 
quasistatic shearing over long total strains produces a 
well defined ensemble of states that is independent of the 
initial configuration, and (ii) the process of shearing, no 
matter how slow, destroys the clustering that can pro- 
duce large </>j's under compression. It is for this reason 

OS 

that (pj represents a true, well-defined, jamming tran- 
sition in the limit of vanishingly small shear strain rate 
7 — >• 0, and does not suffer from the questions of equili- 
bration and protocol that jamming from compression or 
cooling does. 
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IV. DISCUSSION AND CONCLUSIONS 

We can relate our results to two recent experiments. 
Lcchcnault et al. [25 , in experiments on vibrated bi- 
disperse brass disks, interpret their results in terms of 
three relevant densities, <f> g < <fij < <fi a (see their Fig. 1). 
4> g they call the glass/jamming transition where the 
structural relaxation time rapidly grows large on experi- 
mental time scales; we can identify this with the behavior 
in our Fig. [5^,. At <fi a they say that the system reaches the 
fully arrested state] we can identify this as the relatively 
large jamming density one can obtain from slow com- 
pression. For (fig < <fi < (f> a , they say "strong vibration 
can still induce micro-rearrangements through collective 
contact slips that lead to a partial exploration of the 



portion of phase space, restricted to a particular frozen 
structure" and they find a diverging time and length scale 
at a (fij ~ 0.842 within this region; they refer to this <fij 
as the rigidity/jamming transition. We believe this is 
the region where small shear displacements remain pos- 
sible (as illustrated in Fig. [7b at low 7 < 0.05) and that 

— OS 

their <fij corresponds to the <fij of quasistatic shearing. 
In another work by Zhang et al. [26] , a system of disks 
was prepared in a stress-free configuration at a density 
<fi = 0.758, but upon shearing at constant 0, the system 
jammed relatively quickly. The comparatively low value 
of (fi in these experiments, as well as the low average con- 
tact number Z ~ 3 they find at jamming, suggests, as the 
authors say, that friction is playing an important role in 
these experiments. Here we point out, however, that ex- 
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FIG. 8: (color online) (a) Strain averaged cluster strength (716)7, 5 ) ener S v P er particle (E/N)^, and (c) jammed fraction 
(/} 7 , vs total strain 7 at <j> — 0.85 > (pj ; plots (d), (e), and (f) are the same quantities at the lower density (f> = 0.84 < <pj . 
The three curves correspond to simulations starting from three different initial configurations: A is the moderately clustered 
configuration, and B is the highly clustered configuration, of Figs. [6] and [71 while C is a configuration with random initial 
particle positions. The system has N = f024 particles. Note, symbols are displayed on every 2000" 1 data point for = 0.85, 
and every 1000*' 1 data point for cjf> = 0.84, to help with curve identification. 



actly the same behavior may be observed in frictionless 
disks, as illustrated by our Fig. [7] 

To conclude, we have considered various approaches 
to the jamming of 2D disks. Consistent with earlier 
works, our results in Figs. [3^, and [7^, show that a rel- 
atively wide range of jamming densities <f>j are possible 
when compressing mechanically stable configurations or 
when cooling thermally equilibrated configurations. We 
can view compression and cooling as quasi- equilibrium 
processes, since they involve changes in the equilibrium 
variables of cj> and T. We see that, rather than jamming 
being defined at a unique density, the value of <f>j from 
such processes is affected by details of the specific pro- 
tocol, such as compression or cooling rate, the relative 
degree of order (particle clustering) in the initial config- 
urations one starts from, and presumably other details of 
the compression or cooling algorithm. For infinitesimally 
slow compression or cooling, it remains unclear if there 
is a well-defined limiting value of (pj that is lower than 
the <^ max — 0.91 of fully phase separated close packed 
lattices of big and small particles. These observations 
suggest that there is no disordered athermal jamming 



point that is the well defined T — > limit of an equilib- 
rium glass transition T g {<j>) in the (<f>,T) plane. Such an 
equilibrium glass transition, if it exists, must by defini- 
tion be protocol independent, whereas athermal jamming 
via compression or cooling appears not to be. 

However, if one follows Liu and Nagel [J , and moves off 
equilibrium into the phase space of shear driven nonequi- 
librium steady states, then we find that there is a unique 
well-defined athermal jamming transition in the phase 
space (</>, T, 7), located at (cpj , T — 0, 7 — ► 0) (note, it is 
more convenient here to think of the nonequilibrium axis 
as being the shear strain rate 7 rather than the shear 
stress a). We stress that the ((j>,T) plane at 7 —> is 
not the plane of equilibrium; rather it is the plane of 
quasistatically sheared steady states. This is most eas- 
ily seen by noting that for tf> > cpj , the quasitatically 
sheared system has a finite average shear stress a (which 
is just the dynamic yield stress <ty(4>)), whereas in equi- 
librium one would expect to have a = 0. Unlike with 
compression or cooling, the limit of infinitesimally slow 
shearing gives a well defined value cpj clearly less than 
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0max- Unlike compression or cooling, the value of <pj is 
independent of the initial configuration one starts from; 
the process of quasistatic shearing creates a well defined 
ensemble that is independent of the starting configura- 
tion. This is illustrated by the results of our Figs. [7}j8] 
where, even starting from a carefully prepared dense un- 
jammed state with a large degree of particle clustering, 
we find that quasistatic shearing (unlike equilibrium pro- 
cesses) destroys the clustering and, after a finite amount 
of shear, restores one to states typical of the quasistatic 
sheared ensemble. Moreover, we believe that the value 
of <pj is independent of the specific details of the shear- 
ing dynamics, provided one is in the limit of overdamped 

OS 

particle motion; the cpj we report here from quasistatic 
shearing, which involves energy minimization rather than 
a specific particle dynamics, agrees well with the value 
we find from a scaling analysis [57] of shear driven states 
at finite strain rates 7 using Durian's mean-field bubble 
dynamics, Eq. pjl. 

Thus, in our model, the athermal jamming of steady 
state shear driven systems occurs at a nontrivial (i.e. 
9j < '/'max) unique point in the (<f>, T, 7) space that 
is independent of any further details of the shearing pro- 
tocol. Our observation of critical scaling in such shear 
driven flow [T31 [22], leads us to conclude that cpj lo- 
cates a true nonequilibrium critical point in the (</), T, 7) 
phase space. 

The question of whether there can exist a sharp equilib- 
rium glass transition in such simple models remains con- 
troversial. Our results on jamming, reported here, make 
it interesting to speculate that, should such an equilib- 
rium glass transition not exist, there may nevertheless be 
a sharp glass transition when one considers the behavior 
of shear driven steady states. 



the number of such conjugate gradient steps per parti- 
cle, Mft eI /N, for all initial configurations of the RAND 
and QS ensembles, near the ensemble specific jamming 
density. In Fig. [9^i we show results for RAND for a 
system with N = 16384 particles at a packing fraction 
4> = 0.8415 sw <^ AND . In Fig. [9)3 we show results for 



QS for a system with N — 4096 at 



0.8430 



kQS 



We see that the states clearly cluster into two groups: 
those with E/N = I0~ 16 , which correspond to the un- 
jammed states for which our minimization has stopped 
upon reaching our lower cutoff treshhold e cu t , and those 
with larger E/N, corresponding to the jammed states. 
There are exceedingly few states, of negligible statistical 
weight, in the region connecting these two clusters; it is 
these states, with very small but finite energy, which take 
the longest to energy minimize. 




(b) QS 



FIG. 9: Scatter plot of energy per particle E/N vs the num- 
ber of conjugate gradient iteration steps per particle Miter/ N 
needed to achieve minimization at cj> w 4>j. (a) is for RAND 
with N = 16384 particles at = 0.8415; (b) is for QS with 
N = 4096 particles at <f> = 0.8430. The horizontal line at 
E/N — 10 -16 are the unjammed states, where the minimiza- 
tion has stopped upon reaching our cutoff threshold e cu t • 
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Appendix A 

In this section we provide several details concerning 
our simulation methods. We first consider the stopping 
criterion E/N < e cut = 10~ 16 which we use to identify 
unjammed states. Let Mi ter be the total number of con- 
jugate gradient iteration steps needed to achieve energy 
minimization of a particular configuration, according to 
the stopping criteria given in Section II. In Fig. [9] we 
show a scatter plot of the energy per particle E/N vs 



In Fig. 10 we plot the number of conjugate gradient it- 
eration steps per particle M- ltcv /N needed to energy min- 
imize the initial configuragtion, averaged over all initial 
configurations, vs the packing fraction tp, for both RAND 
and QS, for several different system sizes N. Note, since 
each initial state for RAND is completely random, while 
each initial state in QS starts as a small affine shear strain 
from a previously energy minimized state, the number of 
iterations needed for RAND is higher than that needed 
for QS. As N increases, M- ltai {4>,N)/N sharpens up to 
a peak located near <f>j. For RAND, we see that the 
peak value appears to be approaching a constant as N 
increases. For QS, we use A7 = 10~ 4 for sizes N < 2048, 
and for this strain increment, the peak number of iter- 
ations per particle also seems to be approaching a con- 
stant as N increases. For N > 2048, however, we use 
A7 = 10 -5 ; each initial configuration is thus closer to its 
previous energy minimum from the prior strain step, and 
so takes fewer (about a factor of two) iteration steps to 
reach its new energy minimum than for the larger A7. 
While decreasing A7 thus appears to make the energy 
minimization more efficient, as A7 decreases we obvi- 
ously need to run longer to reach the same total strain 
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7- Thus the net effect of decreasing the strain step from 
A7 = 10~ 4 to 10~ 5 is about a factor 5 increase in com- 
putation. It is interesting to note that for RAND, and 
also for QS for those sizes N < 2048 where a constant 
A7 = 1CP 4 is used, the curves of M[ tel (4 l , N) for differ- 
ent N all seem to intersect at roughly the same value 
<f>* k, <pj. If we interpret the total number of iterations 
needed to energy minimize, Mit or , as a relaxation time r, 
then a finite size scaling analysis would suggest a diver- 
gent relaxation time r ~ N ~ L Zag at (j>j, with z cg « 2. 
We note, however, that the conjugate gradient "dynam- 
ics" does not necessarily correspond to a real physical 
dynamics, so z cg need not equal the physical dynamic 
critical exponent z, such as one might find using the dy- 
namics of Eq. pi . 



0.85 




shows that these values are small enough that we are 
in the correct quasistatic limit, except possibly for our 
largest system size N — 8192; we have not been able 
to simulate N = 8192 long enough with A7 < 10~ 5 to 
verify that a smaller A7 is not needed for this case. 
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FIG. 11: (color online) (a) Dependence of jammed fraction 
/ and (b) energy per particle E/N on shear strain step in- 
crement A7, within quasistatic shearing simulations, at fixed 
packing fraction <f> = 0.8430 fts cj)j ■ Results are shown for 
several different system sizes with particle numbers N. 



0.85 



FIG. 10: (color online) Average number of conjugate gra- 
dient iteration steps per particle M- lteT /N needed to energy 
minimize, vs packing fraction <f>. (a) is for RAND; (b) is for 
QS. Results are shown for a few different systems sizes with 
particle numbers N . 



Finally, in Fig. [TT] we consider the effect of the hnite 
shear strain step A7 on our QS simulations. In Fig. [IT]} 
we plot the fraction of jammed states / vs A7 for the 
fixed value <f> = 0.8430 ~ (pj , for systems with different 
numbers of particles N < 4096; in Fig. [Tip we plot the 
corresponding energy per particle E/N vs A7. Note, at 
this value of 0, the jammed fraction / is in the range 
0.5 < / < 0.8, depending on system size N, hence we are 
adequately sampling the behavior in both jammed and 
unjammed states. We see that for fixed N, both / and 
E/N decrease to a constant value (within the estimated 
statistical error |28j ) as A7 decreases. The larger the 
value of N, the smaller A7 must be to reach the limiting 
constant value. For the results reported in the main body 
of this work we have used A7 = 10~ 4 for N < 2048 and 
A7 = 10~ 5 for N > 4096. Comparison with Fig. Ill 
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